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Abstract. In the previous paper an adaptive filtering based on a reference recursive recipe was 
developed and tested on a simulated dynamics of a spring, mass, and damper with a weak nonlinear 
spring. In this paper the above recipe is applied to a more involved case of three sets of airplane data 
which have a larger number of state, measurements, and unknown parameters. Further the flight 
tests cannot always be conducted in an ideal situation of the process noise and the measurement 
noises being white and Gaussian as is generally assumed in the Kalman filter. The measurements 
are not available in general with respect to the center of gravity, possess scale and bias factors 
which will have to be modelled and estimated as well. The coupling between the longitudinal 
and lateral motion brings in added difficulty but makes the problem more interesting. At times 
the noisy measurements from the longitudinal and lateral motion are input into the longitudinal 
states. This leads to the resulting equations becoming linear with the measurement noise forming 
the process noise input. At times it turns out that even a parameter that strongly affects the 
airplane dynamics is estimated which vary widely among the approaches. This requires a careful 
look at the estimates. We also recommend a closer look at the correlation coefficients (that is 
generally ignored in such studies) of the estimated parameters which provide an insight into their 
subsequent uses. The present recipe has been shown to be better than the earlier approaches in 
estimating the unknowns. In particular the generalized cost functions that are introduced in the 
present work help to identify definitive results from deceptive results. 

Keywords. Adaptive EKF, Longitudinal and lateral flight dynamics, Recursive parameter esti¬ 
mation, Cramer Rao Bound. 


1 Introduction 

In the previous paper an extensive study was carried out using an adaptive Extended Kalman 
Filter(EKF) tuning procedure called reference recursive recipe applied to the simulated data of 
a simple spring, mass, and damper system with a weak nonlinear spring. It is useful to further 
demonstrate the effectiveness of the above approach to handle more involved real data studies. It 
turns out that a good choice would be the flight test data analysis of an airplane system dealing 
with longitudinal and lateral motion. These generally have many states and measurements and a 
large number of aerodynamic parameters to be estimated. Further the flight tests cannot always 
be conducted in an ideal situation of the process noise and the measurement noises being white 
and Gaussian as is generally needed in the Kalman filter. The measurements are not available 
with respect to the center of gravity, possess scale and bias factors which will have to be modeled 
and estimated as well. The coupling between the longitudinal and lateral motion brings in added 
difficulty but makes the problem more interesting. At times the noisy measurements from the 
longitudinal and lateral motion are input into the longitudinal states. This leads to the resulting 
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equations becoming linear with the measurement noise forming the process noise input. This is 
another example of subjectivity in estimation theory. However the final results should be mean¬ 
ingful, reasonable, acceptable and useful no matter whatever subjective inputs are introduced into 
the problem formulation and solution. 


2 Aircraft Equations of Motion 

The equations of motion for aircraft flight test data analysis are usually based on the following 
assumptions. The flat non rotating earth forms an inertial system, and with constant gravitational 
acceleration, the atmosphere is fixed relative to earth axes, and the aircraft is rigid with constant 
mass and moment of inertia, and so the effects of fuel sloshing, structural deformations, and the 
relative motion of control surfaces are negligible, further possess a vertical plane of symmetry, 
and the thrust is directed along the longitudinal body axis and through the aircraft center of 
gravity. 


2.1 Aerodynamic Modeling 

The most difficult part of aircraft dynamical analysis is the specification of the aerodynamic forces 
and moments acting on it and consequently most effort in aerospace studies involving analytical, 
computational, wind tunnel and flight test are directed towards it with the last being most realistic 
in terms of aircraft geometry, size and flight conditions. 

The aerodynamic forces and moments depend on the complete history of the aircraft motion in¬ 
volving airspeed, angles of incidence of the air with respect to the aircraft body, its linear and 
angular accelerations, control surface deflections, so a lot of effort goes in determining an ade¬ 
quate and useful aerodynamic model characterization. The aerodynamic modeling is a balance 
between complexity and acceptability for the purpose of control design, simulation, and flying 
qualities. 


2.2 Linearization for System Identification 

Data used for system identification applied to aircraft arise mostly from maneuvers that excite 
the longitudinal short period dynamics or the lateral modes associated with body axis roll, spiral 
motion, and lateral oscillations. In general the flight test data from the longitudinal and lateral 
motions should be decoupled but in practice there could be unavoidable coupling due to the 
airspeed changes and variations in the lateral quantities during a longitudinal maneuver, and vice 
versa. 

For flight data analysis the above coupling can be avoided in two ways. In the first linearization 
is carried out for small motions about a reference condition such as steady, wings level flight with 
no sideslip whence the equations decouple into the longitudinal motion in the plane of symmetry 
and the other the lateral motion out of the plane of symmetry. The resulting linearized models 
are adequate for (i) most flight conditions involving changes in linear and angular velocities from a 
reference condition, (ii) the aerodynamic effects described by linear functions of state and control 
variables, and (iii) the consequent stability analysis and control system design. 

The second method of handling the unavoidable coupling of longitudinal and lateral motion is to 
substitute the measured values in the nonlinear equations in order to linearize them. The resulting 
linear equations are fewer and with less number of model parameters. Here the dynamic pressure 
q , airspeed V, and other variables that cause coupling are substituted by the measured values. For 
flight conditions near stall, spin, and maneuvers involving large angles and angular rates nonlinear 
models must be used. 

2.3 Some Approaches in an Airplane Flight Test Data Analysis 

For airplane flight test data analysis many approaches have been suggested to determine the 
parameters and the noise covariances. These are the natural formulation of Schultz (1976),the 
innovation formulation of Stepner and Mehra (1973), a combined formulation called MMLE3 of 
Maine and Iliff (1981), the combined formulation of Ishimoto (1997) to solve the practical problems 
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in natural and innovation formulations. A more detailed discussion on aircraft flight test data 
analysis can be found in Klein and Morelli (2006) and Jategaonkar (2006). 


2.3a Natural Formulation 

In the natural formulation the innovation cost J is directly minimized using a modified NR tech¬ 
nique with respect to all the unknowns Po, 0, Q, and R. This method requires more computer 
time and memory than an output error method due to Q compulsively bringing in the Kalman 
filter. The computation of all the gradients dependent on the state and measurement functions 
and Q needs more computation. The serious convergence problems in this formulation to estimate 
0 have been illustrated by Maine and Iliff (1981) for a simple scalar problem. 


2.3b Innovation Formulation 

The unknown Q and R appear in J indirectly through innovation covariance and K. The innovation 
formulation trades Q and R for innovation covariance ^ = HPH T + R and K which eliminates the 
convergence and computational problems associated with the estimation of Q. This estimation is 
done as a two step procedure. First an iteration of the NR algorithm is done to revise the estimates 
of all the parameters except for the term HPH T . In the second the estimate of HPH T is revised, 
using the residuals form the previous step and the above steps are repeated until convergence. The 
elements of K are estimated along with the unknowns using the NR algorithm. Being simple is 
most widely used to solve the maximum likelihood estimation problem when both R and Q are 
unknown. 


2.3c Mixed Formulation 

This is a combined formulation, with the best features of the natural and innovation formulations 
in minimizing J. Instead of Po, Q, and R one works with innovation covariance 2%, K, and R but 
should ensure a consistent K which can be in general be not square like innovation covariance and 

R. 

2.3d Combined Formulation 

This solves the practical problems of natural and innovation formulations and different from 
MMLE3. The unknown parameters are assumed to be included in the system transition matrices 
as well as Q and R. The filter provided innovation covariance 2%f is compared with the ones com¬ 
puted from the residuals 2% r . In this method J is minimized subject to C(0) = diag(^/ — 2& r )=0. 
The purpose of this constraint is to accelerate the estimation of R. Although the sample covari¬ 
ance 21f is a function of the unknown parameters it is fixed during each optimization step and it is 
updated after every iteration. The Sequential Quadratic Programming method used to solve the 
above optimization problem provided substantial improvement in convergence. 

One may ask the question as to why there are so many formulations for solving an optimization 
problem. The reason is the unknowns do not occur in a simple way in the cost function, and there 
are many transformed variables with which one tries to solve for the basic unknowns. Further the 
size and the required compatibility conditions among the transformed variables lead to the many 
difficulties not found in the classical optimization problems. 


3 Analysis of Real Flight Test Case - 1 

The salient features of this aircraft are available in NASA TM-X 56036 (Shafer 1975) and in NASA 
TP 1690 (Maine and Iliff 1981). The parameters are estimated in dimensionless form. The data 
set obtained is for a short period motion excited by the up and down elevator control input ( S e in 
degrees) as shown in Fig.l. In general the flight test data is such that the longitudinal and lateral 
motions are decoupled. Some of the available measurements have been used as inputs in the state 
equations which includes roll angle (</> m ), sideslip (/3 m ), roll rate (p m ), yaw rate (r m ) and the angle 
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of attack (a m ) are shown in Fig.2, Fig.3, Fig.4, Fig.5 and Fig.9 respectively. The state equations 
(n = 3) for the angle of attack (a), pitch rate (q) and the pitch angle (0) respectively are 


—£1 

a =- —C L + q + fz{cos{4> m )cos{a m )cos(6) + sin{a m )sin{6))- 

mV V 

/3 m (p m cos(a m ) + r m sin(a m )) 

qSc c Izz — Ixx 

q — . \Cm a T Cm q TZTpq T C ms 0 e T 6 mo ) - VmPm 

Iyy 2V ' Iyy 

9 =qcos(4> m ) - rmsin^m) + 9 0 
The measurement equations (m=5) are 

o; m — Of ^ 

Qm = Q 
dm = 9 

_ qS n , Xa n . 

a n m — - CN 4- q 

mg g 

qS z ax . 

a Xm — - Ca H- q 

mg g 


where 


Cl =Cncos{q) — CAsin(a) + Cl 0 

Cn = Civ Q + Cjv 6e + Cn 0 

C'A =CA a OL + CA a 2 + CA s ^S e + C^o 

The unknown parameter set (p = 13) is 0 = ( C Na , C Nde , C Lo , C m „, C mq , C mSe , C mo , 9 0 , C No , C Aa , 
C’A 2 > Ca s , Ca 0 ) T ■ The ones with suffix l (5 e ’ are the control derivatives, the ones with suffix zero are 
the biases and all others are aerodynamic derivatives. The initial states are taken as initial measure¬ 
ment and the initial parameter values are taken as (4, 0.24, 0.17, —0.48, —17, —0.9, —0.05, —0.02, 
0.175, -0.3,0.03, -0.083, -0.015) T . 

Other constant values used for case-1 


c=5.58 

00 

II 

CO 

m=172.667 

Ixx=4142.9 

Iyy=3922.4 

Izz=7642.5 

g=32.2 

V=403.1 

<7=83.08 

K a x a =-0.0279 

®a„=0.101 

*».=-1-17 


In the analysis of real flight data the following filter outputs are studied which provide an insight 
into the filter performance. The convergence of the following quantities through the iterations are 
analyzed. 

1. The parameter estimates 0 and their covariances P©. 

2. The noise covariance of Q and R. 

3. The state dynamics without measurement and process noises based on the estimated param¬ 
eter after the filter pass through the data Xd, the prior state X-, the posterior state X+, the 
smoothed state Xs and the measurement Z. 

4. The sample innovation, filtered residue and the smoothed residue along with their bounds 
which is the square root of the predicted covariances given respectively by (R +Hk.Pk\k-\H'^), 
(R - Pfc|fcPfc|fcPfj fc ) and (R - H k \ N P k \ N H^ N ) by the filter. 

5. The estimated measurement and process noise samples as well as their autocorrelations. 

6. The cost functions (J1-J8) after the final convergence. 

Case-1 real data is run using the reference EKF (Q > 0) with 100 iterations. It turns out that 
the off diagonal elements of the correlation coefficient matrix with Q = 0 reduced substantially 
when Q > 0 in all the present case studies thus indicating the latter estimates more trustworthy. 
The Fig.1-5 show the inputs used in state equations. The Fig.6 shows the variation of parameter 
estimates and its initial covariance Po with iterations and a similar Fig.7 for Q and R. The values 
of J1-J3 are close to the number of measurements (to = 5) with J6-J8 are close to the number of 
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states (n = 3) as shown in Fig.8 and Table-1. This means the measurement and state equations are 
well balanced. The J5 is the negative log likelihood cost function. The later Fig.9-10 compares (i) 
the state dynamics based on the estimated parameter after the filter pass through the data, (ii) the 
state after measurement update, (iii) the smoothed state and (iv) the measurement. Unlike in the 
simulated studies, the estimated measurement and process noise did not have constant statistical 
characteristics across time. Another experiment was carried out by generating a typical data set 
by using the estimated theta and injecting the estimated Q and R as additive white Gaussian 
noise. This is to determine the effect of non White and non Gaussian noise distribution in the 
real data on the CRBs. After each iteration in the reference recipe the 0, Q and R were reset 
as from the real data. Similar experiment was also conducted by updating 0 as well. It turned 
out that there is not much of a difference in the final estimates and the CRBs. Two other filter 
runs were carried out using the MT and MS statistics for the estimation of Q and R with scaled 
up Po- The behaviour of the various cost function and in particular J6 and J7 in Table-2 shows 
that the choice of the filter statistics for estimating Q and R in the proposed reference approach 
is the best possible when compared to other approaches. Another feature of recursive parameter 
estimation is that it can vary through time instants and point to two distinct values as reflected 
in the estimation of Cjsr a in Fig.11. This feature of tracking time varying parameters by the EKF 
brings in clearly another advantage of sequential processing instead of batch processing of the data 
by a least squares (LS) procedure. If LS had been used then for the parameter ‘CnJ only an 
average value would have been obtained. If one persists in using a batch processing procedure 
then to get the varying parameter such a feature should be modelled which would include where 
to change and one might have to use the data in blocks by splitting them and all such exercises 
have to be carried out which are not easy. 


3.1 Remarks on Case - 1 

The NASA results have been generated assuming Q = 0 and are comparable with reference pro¬ 
cedure for the parameter estimates and their CRBs. Further the MT and MS methods give quite 
different estimates for the Q and R values than from the reference case. We believe that the 
reference procedure provides the best possible parameter estimates and their uncertainties. The 
percentage CRB of the i th parameter estimate is defined as %CRB=100 x y / Pe../|0 ?; |. From the 
plot of % error in the parameter estimate with respect to the reference value and the %CRB in 
Fig 12, it can be seen that the parameters f? 0 and Ca & are relatively weak when compared to 
other parameters. The CRBs as estimated by different methods generally appear to vary widely. 
However what is interesting is that even the estimate of the strong parameter C mq varies widely 
among the methods. Such a behaviour of the filter across the parameter estimates shows how im¬ 
portant is the tuning of the filter statistics namely Po, Q and R in parameter estimation and their 
uncertainties. The rounded off 100xC matrix of the parameter estimates for this case is 


0 

C Na 

Cn Sf 

Cl 0 

Cm a 

Cm q 

Cms r 

Cm 0 

o 0 

C N 0 

C Aa 


C ASf 


c Na 

100 

62 

0 

-18 

0 

-7 

-3 

0 

32 

-1 

-10 

-10 

-8 

c NS( 

62 

100 

0 

-10 

1 

-11 

-10 

0 

93 

-4 

-2 

-16 

-15 

Cl 0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Cm a 

-18 

-10 

0 

100 

18 

51 

31 

0 

-5 

-15 

1 

-15 

-7 

Cm q 

0 

1 

0 

18 

100 

76 

79 

0 

1 

6 

-3 

5 

2 

n 

-7 

-11 

0 

51 

76 

100 

97 

0 

-10 

-1 

-2 

-14 

-14 

C mo 

-3 

-10 

0 

31 

79 

97 
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0 

-11 

2 

-3 

-11 

-15 

0o 

0 
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0 

0 
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0 

0 

0 

CjVo 

32 

93 
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Figure 1: Control input (5 e in degrees) versus time (s) 



Figure 2: Measurement input {(j> m in degrees) versus time (s) 



Figure 3: Measurement input (/3 m in degrees) versus time (s) 



Figure 4: Measurement input ( p m in degrees) versus time (s) 



Figure 5: Measurement input (r m in degrees) versus time (s) 
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Figure 9: Comparison of the predicted dynamics, posterior, smoothed and 
measurement in degrees (a. angle of attack b. pitch rate c. pitch angle) vs time. 
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Figure 10: Comparison of the predicted dynamics, posterior, smoothed and 
measurement in ft/sec 2 (a.normal acceleration b.axial acceleration) vs time 



Figure 11: Variation of the parameter estimate (Cjy a ) through time instants 
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Figure 12: Comparison of (a) Absolute percentage error with respect to the reference of the 
parameter estimates and (b) %CRBs by different methods 


Table 1: Real flight test data case-1 results (0,<re)- 


0 

Reference 

NASA 

Gemson 

MT 

MS 

c Na 

4.6469 

4.9584 

4.7073 

4.6978 

4.9141 

(0.0179) 

(0.1168) 

(0.039) 

(0.0229) 

(0.0422) 

C Nse 

0.0555 

0.3023 

0.1292 

0.1225 

0.4691 

(0.0277) 

(0.1550) 

(0.0523) 

(0.0357) 

(0.0517) 

Cl 0 

0.0162 

0.2189 

-0.0064 

0.0160 

0.0184 

(0.0032) 

(0.009344) 

(0.0048) 

(0.0018) 

(0.0021) 


-0.5468 

-0.6125 

-0.63 

-0.5560 

-0.5885 

(0.0093) 

(0.00953) 

(0.0188) 

(0.0098) 

(0.0036) 

Cm q 

-19.8027 

-22.27 

-20.8623 

-19.7062 

-20.2395 

(0.6692) 

(0.7713) 

(1.1908) 

(0.7286) 

(0.2937) 

n 

ms e 

-1.1229 

-1.2193 

-1.2763 

-1.1396 

-1.1503 

(0.0218) 

(0.02881) 

(0.0442) 

(0.0236) 

(0.0111) 

Cm 0 

-0.0495 

-0.0532 

-0.0561 

-0.0502 

-0.0497 

(0.0012) 

(0.00165) 

(0.0023) 

(0.0013) 

(0.0006) 

r\ 

0.0007 

0.0273 

0.0007 

0.0008 

0.0003 

Vo 

(0.0021) 

(0.04518) 

(0.0135) 

(0.0011) 

(0.0012) 

Cn 0 

0.2195 

0.2254 

0.2225 

0.2218 

0.2358 

(0.0014) 

(0.008725) 

(0.0029) 

(0.0018) 

(0.0028) 

C Aa 

-0.1398 

-0.3639 

0.0214 

-0.1401 

-0.1265 

(0.0153) 

(0.05328) 

(-0.1023) 

(0.0185) 

(0.0197) 

C A „ 2 

-3.2088 

- 

-3.2397 

-3.2088 

-3.8625 

(0.1702) 

(-) 

(0.2430) 

(0.2070) 

(0.2376) 

C Ase 

-0.0651 

-0.07 

-0.0267 

-0.0633 

-0.1178 

(0.0134) 

(0.08084) 

(0.0191) 

(0.0160) 

(0.0167) 

Ca 0 

-0.0155 

-0.0131 

-0.0144 

-0.0154 

-0.0182 

(0.0007) 

(0.004088) 

(0.0010) 

(0.0008) 

(0.0008) 


Table 2: Real flight test data case-1 results* (R,Q,J). 


R (Ref) 
x 10 -6 

Q (Ref) 
xlO -6 

J1-J8 

(Ref) 

R (MT) 
x 10 -6 

Q (MT) 
x 10 -6 

J1-J8 

(MT) 

R (MS) 
x 10 -6 

Q (MS) 
x 10 -6 

J1-J8 

(MS) 

0.49 

0.04 

0.40 

15.98 

17.70 

0.134 

2.287 

1.204 

4.4752 

5.1532 

4.6432 

0.0004 

-56.2206 

2.9551 

2.9303 

2.5161 

0.4107 

0.0312 

3.9381 

94.5086 

26.3511 

0.0393 

2.6418 

0.3231 

4.0090 

3.9630 

2.9764 

0.0004 

-54.7596 

6.6681 

6.4985 

2.4562 

3.2046 

37.6770 

7.5509 

198.2716 

28.9841 

0.0001 

0.0015 

0.3456 

3.3893 

3.3866 

3.2057 

0.0002 

-49.6223 

3.8921 

4.7110 

2.6369 


*Cost functions are not close to their expected values in MT and MS methods. 
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4 Analysis of Real Flight Test Case - 2 


The data set is obtained from NASA TP 1690 (Maine 1981) by employing a peculiar manoeuvre 
where elevator control input (S e in degrees) shown in Fig. 13 is imparted when the aircraft (T 37 B) 
is rolling through a full rotation about its x-axis during aileron roll. Similar to the earlier case, the 
coupling between the longitudinal and lateral motion are replaced by their measured values which 
includes roll angle ( <f > m ), sideslip (/3 m ), velocity ( V m ), roll rate ( p m ), yaw rate (r m ) and the angle 
of attack ( a m ) as shown in Fig.15, Fig.16, Fig.17, Fig.18, Fig.14 and Fig.22 respectively. The state 
equations (n = 3) for the angle of attack (a), pitch rate (q) and the pitch angle (9) respectively 
are 


a = — 


qS 


' ol + Cl s S e + Cl 0 ) + q + 


9 


-(cos(4> rn )cos(a m )cos(6)+ 


mV m cos(/3 m ) “ Se 0 V m cos{!3 n 

sin(a m )sin(6 )) - tan(/3 m )(p m cos(a m ) + r m sin(a m )) 

qSc,„ ^ c c . Izz — Ixx 

9 — T (G m & T Cm. XTT 9 lC mi XTyd C'm; O e + C'm.Q ) T 7 X m Pm 
lyy q 2V 2V e lyy 

9 =qcos{(f> m ) - r m sin(4> m ) + 9 0 
The measurement equations (m=4) are given by 


‘ V 

9m = 9 

9m = 9 

qQ X 

a n m = {Cjsi a + Cn Sb S e + Cn 0 ) 4- —q 

mg g 


The unknown parameters (p = 10) is ( C La , C Lse , C Lo , C m<x , C mq , C m& , C m5e , C mo , 9 0 , C No ) T with 
an approximation C^ a = C^ a and Cjv 4 = Cr J/s . The ones with suffix ‘<5 e ’ are the control 
derivatives, the ones with suffix zero are the biases and all others are aerodynamic derivatives. 
The initial states are taken as initial measurement and the initial parameter values are taken as 
(4,0.15,0.2, -0.5, -11.5, -5, -1.38, -0.06, -0.01,0.2) T . 


Other constant values used for case-2 


S=184 

m=196 

Ixx=6892.7 

Iyy=3953.2 

Izz=10416.4 

g=32.2 

c=5.58 

K a x a =- 0.0279 

®a„=0-101 

I< a = 1 


Case-2 real data is run using the reference EKF (Q > 0) with 100 iterations. The Figs.13-18 are 
the inputs used in state equations. The Fig. 19 shows the variation of parameter estimates and 
its initial covariance Po with iterations and a similar Fig.20 for Q and R. The values of J1-J3 
are close to the number of measurements (to = 4) with J6-J8 are close to the number of states 
(n = 3) as shown in Fig.21 and Table-3. This means the measurement and state equations are well 
balanced. The J5 is the negative log likelihood cost function. The later Figs.22-23 compares (i) 
the state dynamics based on the estimated parameter after the filter pass through the data, (ii) the 
state after measurement update, (iii) the smoothed state and (iv) the measurement. Unlike in the 
simulated studies, the estimated measurement and process noise did not have constant statistical 
characteristics across time Another experiment was carried out by generating a typical data set 
by using the estimated theta and injecting the estimated Q and R as additive white Gaussian 
noise. This is to determine the effect of non white and non Gaussian noise distribution in the real 
data on the CRBs. After each iteration in the reference recipe the 0, Q and R were reset as from 
the real data. Similar experiment was also conducted by updating 0 as well. It turned out that 
there is not much of a difference in the final estimates and the CRBs. Two other filter runs were 
carried out using the MT and MS statistics for the estimation of Q and R with scaled up Po- 
The behaviour of the various cost function and in particular J6 and J7 in Table-4 shows that the 
choice of the filter statistics for estimating Q and R in the proposed reference approach is the best 
possible when compared to other approaches. 
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4.1 Remarks on Case - 2 


The NASA results have been generated assuming Q = 0 and are comparable with reference pro¬ 
cedure for the parameter estimates and their CRBs. Further the MT and MS methods give quite 
different estimates for the Q and R values than from the reference case. We believe that the 
reference procedure provides the best possible parameter estimates and their uncertainties. From 
the plot of % error in the parameter estimate with respect to the reference value and the %CRB 
in Fig 24, it can be seen that the parameters Cl s and are relatively weak when compared to 
other parameters. The CRBs as estimated by different methods generally appear to vary widely. 
However what is interesting is that even the estimate of the strong parameter such as C mq varies 
widely among the methods. Such a behaviour of the filter across the parameter estimates shows 
how important is the tuning of the filter statistics namely Po, Q and R in parameter estimation 
and their uncertainties. It was also observed that for this particular case, the cost J2 using filtered 
residue went negative at some iteration as seen in Fig.21 with a spike whose absolute value was 
plotted. However the cost J3 using smoothed residue that was used for tuning the filter did not 
show any such peculiarity. The rounded off 100 xC matrix of the parameter estimates for this case 
is 
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Figure 13: Control input (6 e in degrees) versus time (s) 



Figure 14: Measurement input (<j) m in degrees) versus time (s) 
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Figure 22: Comparison of the predicted dynamics, posterior, smoothed 
and the measurement in degrees(a. angle of attack b. pitch rate c. pitch angle) vs time 
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Figure 23: Comparison of the predicted dynamics, posterior, smoothed 
and the normal acceleration measurement in ft/seer vs time 
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Figure 24: Comparison of (a) Absolute percentage error with respect to the reference of the 
parameter estimates and (b) %CRBs by different methods 


Table 3: Real flight test data case-2 results (0,<re)- 


0 

Reference 

NASA 

Gemson 

MT 

MS 

c La 

4.9235 

5.1068 

4.9028 

4.9260 

5.0620 

(0.0164) 

(0.1322) 

(0.0168) 

(0.0184) 

(0.0323) 

C Lse 

0.1554 

0.1909 

0.0879 

0.1587 

0.3594 

(0.0271) 

(0.1602) 

(0.0267) 

(0.0302) 

(0.0508) 

Cl 0 

0.2409 

0.2448 

0.2529 

0.2408 

0.2517 

(0.0021) 

(0.009215) 

(0.0018) 

(0.0023) 

(0.0027) 

Cm* 

-0.5293 

-0.6474 

-0.6174 

-0.5285 

-0.5590 

(0.0079) 

(0.02339) 

(0.0211) 

(0.0082) 

(0.0055) 

Cm q 

-11.8596 

-14.26 

-18.8339 

-11.8255 

-12.5965 

(0.2402) 

(0.6528) 

(0.8379) 

(0.2483) 

(0.1400) 

C m& 

-6.8959 

-8.27 

-7.1290 

-6.8798 

-6.6713 

(0.4891) 

(1.296) 

(1.544) 

(0.5062) 

(0.3021) 

r 

y ~ / rns e 

-0.9731 

-1.1614 

-1.1841 

-0.9711 

-1.0247 

(0.0177) 

(0.05371) 

(0.471) 

(0.0184) 

(0.0129) 

Cm 0 

-0.0425 

-0.0505 

-0.0507 

-0.0424 

-0.0447 

(0.0009) 

(0.002655) 

(0.0024) 

(0.0009) 

(0.0006) 

A 

0.0003 

-0.01177 

-0.0037 

0.0002 

-0.0006 

C'O 

(0.0021) 

(0.02528) 

(0.001) 

(0.0011) 

(0.0007) 

Cn 0 

0.2538 

0.2541 

0.2503 

0.2540 

0.2635 

(0.0014) 

(0.008935) 

(0.0014) 

(0.0016) 

(0.0026) 


Table 4: Real flight test data case-2 results* (R,Q,J). 


R (Ref) 
xlO" 6 

Q (Ref) 
xl0~ 6 

J1-J8 

(Ref) 

R (MT) 
xlO" 6 

Q (MT) 
xlO" 6 

J1-J8 

(MT) 

R (MS) 
xlO" 6 

Q (MS) 
x 10“ 6 

J1-J8 

(MS) 

1.241 

0.051 

0.460 

5.668 

0.180 

2.954 

2.646 

3.9336 

4.2225 

3.6162 

0.0008 

-44.1347 

2.9752 

2.9760 

2.9070 

1.6135 

0.2395 

2.3155 

2.9290 

0.2025 

3.1532 

0.6666 

3.7662 

4.5191 

3.8384 

0.0008 

-43.7340 

4.2266 

4.2284 

2.9489 

3.1599 

37.2424 

9.3413 

841.5496 

0.00005 

0.0003 

0.2386 

3.1621 

3.1507 

2.5900 

0.0007 

-38.0517 

8.4768 

8.4655 

3.0215 


*Cost functions are not close to their expected values in MT and MS methods. 
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5 Analysis of Real Flight Test Case - 3 


The data set is obtained from NASA TP 1690 which describes the lateral motion of a oblique wing 
aircraft with zero wing skew excited by the control input (<5 0 and S r in degrees) as shown in Fig.25. 
Similar to the earlier case, the coupling between the longitudinal and lateral motion are replaced 
by their measured values which includes pitch angle (9 rn in rad), pitch rate (q rn in rad/s) and the 
angle of attack (a m in rad) as shown in Fig.26, Fig.27 and Fig.28 respectively. The state equations 
(n=4) for the angle of sideslip (/3), roll rate (p), roll angle (4) and yaw rate (r) are 


P — ~y (Cy^/3 + C Yp ^yP + Cy T ^yr + C Yfa S a + C Ysr S r + Po)+ 
y sin(4)cos{9 m ) + psin{a m ) - rcos(a m ) 

P ~ r TZZ + Cl p -jyP + Cl p -^yT + CL Sa S a + C Ltr S r + Cl 0 )+ 


IXX IXX 

Iyy — Izz 


Izx 


T -rq m + t— pq m 

1XX 1XX 


4 =p + q tan(9 m )sin((i i) + r tan(9 m )cos(4) + 4o 

Ixx — Iyy 


t — F C N p~2yP F ^ Nr 2Y r F ^N 5a S a + CN 5r S r + Cn 0 )+ 


Izx 

T -PQm - 7— rq m 

Izz Izz 


The measurement equations (m= 5) are given by 


n ' P p r 
fim Zy — T KyXy — 

Pm =P 

= 4* 


r m = r 

a Vm = + Cy p —P + ( -'Y r -^y r + CYs a S a + 0 Ysr S r + Cy 0 )- -^~P ~\ - ~^-T 


The unknown parameter set (p = 20) is 0 = (C Yf) , C Y g r , Po, C Lft , C Lp , C Lr , C LSa , C L , r , C Lo , 4 o, C Ng , 
Cn p , Cn p , Cpf Sa , Cn Sp , Cn 0 , C Yo , Cy- p , C Yr , C Ys ) T ■ The ones with suffix ‘(5 a ’ and ‘S r ' are the con¬ 
trol derivatives, the ones with suffix zero are the biases and all others are aerodynamic derivatives. 
The initial states are taken as initial measurement and the initial parameter values are taken as 
(-0.5,0.1, -0.01,0.01, -0.35,0.01,0.06,0.01, -0.002,0.002,0.07, -0.055, -0.05,0.003, -0.04,0.0068, 
- 0.025,0.5,-1,0.005) t . 


Other constant values used for case-3 


q = 865.3 

S=9.3 

m=387.7 

Ixx=314 

Iyy=488 

Izz=698 

Izx=69 

V=39.41 

g=9.81 

b=6.81 

K p zp= 0.305 

Kyxp=2.7?> 

Ai„ =-0-098 

x av =0.651 


Case-3 real data is run using the reference EKF (Q > 0) with 100 iterations. Fig.25-28 are the 
inputs used in state equations. The Fig.29 shows the variation of parameter estimates and its 
initial covariance Po with iterations and a similar Fig.30 for Q and R. The values of J1-J3 are 
close to the number of measurements (to = 5) with J6-J8 are close to the number of states (n = 4) 
as shown in Fig.31 and Table-5. The J5 is the negative log likelihood cost function. The later 
Fig.32-34 compares (i) the state dynamics based on the estimated parameter after the filter pass 
through the data, (ii) the state after measurement update, (iii) the smoothed state and (iv) the 
measurement. Unlike in the simulated studies the estimated measurement and process noise did 
not have constant statistical characteristics across time. Another experiment was carried out by 
generating a typical data set by using the estimated theta and injecting the estimated Q and R 
as additive white Gaussian noise. This is to determine the effect of non white and non Gaussian 
noise distribution in the real data on the CRBs. After each iteration in the reference recipe the 0, 
Q and R were reset as from the real data. Similar experiment was also conducted by updating 0 
as well. It turned out that there is not much of a difference in the final estimates and the CRBs. 
Two other filter runs were carried out using the MT and MS statistics for the estimation of Q 
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and R with scaled up Po- The behaviour of the various cost function and in particular J6 and 
37 in Table-6 shows that the choice of the filter statistics for estimating Q and R in the proposed 
reference approach is the best possible when compared to other approaches. 


5.1 Remarks on Case - 3 

The NASA results have been generated assuming Q = 0 and are comparable with reference pro¬ 
cedure for the parameter estimates and their CRBs. Further the MT and MS methods give quite 
different estimates for the Q and R values than in the reference case. We believe that the reference 
procedure provides the best possible parameter estimates and their uncertainties. From the plot of 
% error in the parameter estimate with respect to the reference value and the %CRB in Fig 35, it 
can be seen that the parameters /3q, <j> o, C N Sa , CV and Cy r are relatively weak when compared to 
other parameters. The CRBs as estimated by different methods generally appear to vary widely. 
However what is interesting is that even the estimate of the strong parameter such as Cl v varies 
widely among the methods. Such a behaviour of the filter across the parameter estimates shows 
how important is the tuning of the filter statistics namely Po, Q and R in parameter estimation 
and their uncertainties. The rounded off 100xC matrix of the parameter estimates for this case 
is 
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In analyzing the correlation coefficient matrix it is firstly useful to see the combination of parame¬ 
ters that occur in the governing equations of the flight data analysis. The sets in the present case 
are the trim values (/3o and $o) have very little correlation with all other estimated parameters. 
The parameters in the three sets (C Lfj , C Lp , C Lr , C Lga , C L/ir , C Lo ), {C Nf> , C Np , C Nr , C Nsa , C Nsr , 
Cn 0 ), and (Cy,, CV , Cy r , Cy , Cy s , Cy 0 ), all have very similar correlations among themselves 
as seen in the blocks of matrices and the reason is as follows. There is coupling of the dynamical 
motion due to the states and the controls. If a certain state or control is excited relatively higher 
than others then the estimated parameter that multiplies it will have lower correlation with other 
estimates in the set and vice versa. Since the parameter sets are similarly excited all of them have 
similar correlation coefficient matrices. This feature is similar to the spring, mass, and damping 
system considered in Part-1 of the paper. 
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Figure 29: Variation of initial parameters Oq (continuous) and its P 0 (dashed) with iterations 
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Figure 30: Variation of Q (dashed) and R(continuous) with iterations 
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Figure 31: Variation of different costs (J1-J8) with iterations 
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Figure 32: Comparison of the predicted dynamics, posterior, smoothed 
and the measurement in degrees (a. Sideslip b. Roll angle) vs time 



Figure 33: Comparison of the predicted dynamics, posterior, smoothed 
and the measurement in degrees/sec (a. Roll rate b. Yaw rate) vs time 



Figure 34: Comparison of the predicted dynamics, posterior, smoothed 
and the measurement in ft/sec 2 (lateral acceleration) vs time 
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Figure 35: Comparison of (a) Absolute percentage error with respect to the reference of the 
parameter estimates and (b) %CRBs by different methods 


6 Conclusions 

The present adaptive filtering approach based on the reference recursive recipe is applied to the 
more involved cases of three sets of real airplane flight test data which have a larger number of 
state, measurements, and unknown parameters. A closer look at the correlation coefficients that 
is generally ignored in such studies of estimating the unknown parameters indicates the necessity 
to process the data by including the process noise Q in addition to the measurement noise R. 
Generally the parameter estimates across the various approaches are close but their CRB can vary 
much more among them. In particular the generalized cost functions based on balancing the state 
and measurement equations using the many filter outputs introduced in the present work help 
to show the present approach to be better than the earlier approaches and also help to identify 
definitive results from deceptive results. 
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Table 5: Real flight test data case-3 results (0, <r©). 


0 

Reference 

NASA 

Gernson 

MT 

MS 

C Y, 

-0.4579 

-0.4792 

-0.4761 

-0.4541 

-0.4642 

(0.0043) 

(0.01711) 

(0.0043) 

(0.0053) 

(0.0049) 

Cy s „ 

0.1040 

0.0887 

0.0981 

0.0741 

0.0797 

(0.0067) 

(0.01955) 

(0.0065) 

(0.0065) 

(0.0057) 

Ad 

-0.0143 

-0.10116 

-0.0124 

-0.0107 

-0.0109 

(0.0048) 

(0.00294) 

(0.0021) 

(0.0034) 

(0.0009) 

Cl, 

-0.0168 

-0.0205 

-0.0182 

-0.0170 

-0.0177 

(0.0005) 

(0.00107) 

(0.0011) 

(0.0004) 

(0.0003) 

Cl p 

-0.3100 

-0.36 

-0.3585 

-0.3112 

-0.3080 

(0.0028) 

(0.00713) 

(0.0048) 

(0.0027) 

(0.0022) 

c Lr 

0.0740 

0.0697 

0.0731 

0.0733 

0.0757 

(0.0030) 

(0.005884) 

(0.0066) 

(0.0028) 

(0.0022) 

c LSa 

0.0557 

0.0612 

0.0622 

0.0557 

0.0546 

(0.0004) 

(0.001050) 

(0.0007) 

(0.0004) 

(0.0003) 

c Lsr 

0.0072 

0.006 

0.0089 

0.0073 

0.0082 

(0.0007) 

(0.001252) 

(0.0031) 

(0.0007) 

(0.0005) 

Cl o 

-0.0020 

-0.002 

-0.0023 

-0.0020 

-0.0021 

(0.0001) 

(0.0001467) 

(0.0003) 

(0.0001) 

(0.0001) 

4> o 

0.0018 

0.1506 

0.0023 

0.0019 

0.0019 

(0.0027) 

(0.07034) 

(0.0011) 

(0.0012) 

(0.0013) 

Cm, 

0.0656 

0.0705 

0.0703 

0.0657 

0.0662 

(0.0005) 

(0.000478) 

(0.0009) 

(0.0004) 

(0.0004) 

Cn p 

-0.0429 

-0.046 

-0.0557 

-0.0473 

-0.0596 

(0.0031) 

(0.004006) 

(0.0039) 

(0.0023) 

(0.0021) 

c Nr 

-0.0880 

-0.1062 

-0.0576 

-0.0854 

-0.1021 

(0.0033) 

(0.003562) 

(0.0045) 

(0.0023) 

(0.0010) 

°N Sa 

0.0004 

0.0006 

0.0033 

0.0010 

0.0030 

(0.0005) 

(0.0005924) 

(0.0006) 

(0.0003) 

(0.0003) 

c Nsr 

-0.0478 

-0.0513 

-0.048 

-0.0476 

-0.0502 

(0.0008) 

(0.0009139) 

(0.0013) 

(0.0006) 

(0.0003) 

Cn 0 

0.0067 

0.0072 

0.0068 

0.0067 

0.0069 

(0.0001) 

(0.0001181) 

(0.0002) 

(0.0001) 

(0.00003) 

Cy 0 

-0.0259 

-0.0242 

-0.0251 

-0.0221 

-0.0228 

(0.0008) 

(0.002307) 

(0.0007) 

(0.0008) 

(0.0007) 

C'Yp 

-0.2828 

- 

- 

0.0821 

0.5223 

(0.0327) 

- 

- 

(0.2941) 

(0.0710) 

Cy, 

0.2224 

- 

- 

-1.2336 

-0.8452 

(0.0281) 

- 

- 

(0.4470) 

(0.1355) 

Cys q 

0.0384 

- 

- 

0.0011 

-0.0017 

(0.0047) 

- 

- 

(0.0026) 

(0.0024) 


Table 6: Real flight test data case-3 results* (R,Q,J). 


R 

xlO” 6 

(Ref) 

Q 

xlO" 6 

(Ref) 

J1-J8 

(Ref) 

R (MT) 
xlO -6 

Q (MT) 
x 10~ 6 

J1-J8 

(MT) 

R (MS) 
xlO" 6 

Q (MS) 
xKT 6 

J1-J8 

(MS) 

0.0871 

0.0623 

0.2255 

0.0200 

43.8064 

4.2163 

5.1340 

4.9426 

1.4324 

4.7650 

4.8321 

3.5272 

0.0004 

-55.0111 

3.9673 

3.9669 

3.8171 

2.83 

18.86 

3.88 

4.09 

73.6 

2.0481 

3.7876 

1.0057 

0.5502 

4.3450 

4.3888 

3.1039 

0.0003 

-51.3490 

9.3006 

9.3005 

3.5105 

13.03 

88.07 

4.5 

36.36 

60.13 

0.0005 

0.0007 

1.0975 

0.0016 

4.8127 

4.8200 

4.5173 

0.0003 

-47.2441 

7.5931 

7.5896 

3.9681 


*Cost functions are not close to their expected values in MT and MS methods. 
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